#include <StVenKir.hpp>

namespace Life
{

  template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
  void
  StVenantKirchhoff<Dim, Order, Entity>::init()
  {
   M_bdf_displ_struct->shiftRight( *U_displ_struct );
   M_bdf_velocity_struct->shiftRight( *U_velocity_struct );
   M_bdf_acceleration_struct->shiftRight( *U_acceleration_struct );

   M_bdf_displ_struct->initialize( *U_displ_struct );
   M_bdf_velocity_struct->initialize( *U_velocity_struct );
   M_bdf_acceleration_struct->initialize( *U_acceleration_struct );

   M_bdf_displ_struct->start();
   M_bdf_velocity_struct->start();
   M_bdf_acceleration_struct->start();
  }

  template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
  void
  StVenantKirchhoff<Dim, Order, Entity>::updateBdf()
  {

      // save data
      this->application()->changeRepository( boost::format( "%1%/%2%/P%3%/h_%4%/" )
                                             % this->application()->about().appName()
                                             % convex_type::name()
                                             % Order
                                             % this->application()->vm()["hsize"].template as<double>()
                                             );
#if 0
      M_bdf_displ_struct->shiftRight( *M_identity_ale );
      M_bdf_velocity_struct->shiftRight( M_aleFactory->getDisplacement() );
      M_bdf_acceleration_struct->shiftRight( *U_fluid );

      M_bdf_displ_struct->next();
      M_bdf_displ_struct->next();
      M_bdf_acceleration_struct->next();
#endif
}


  template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
  void
  StVenantKirchhoff<Dim, Order, Entity>::run()
  {
    using namespace Life::vf;

    //----------------------------------------------------------------------------------------------------------//

    if ( this->application()->vm().count( "help" ) )
      {
	std::cout << this->application()->optionsDescription() << "\n";
	return;
      }

    //----------------------------------------------------------------------------------------------------------//

    // First we create the mesh
    //mesh_ptrtype mesh = createMesh( meshSize );

    //----------------------------------------------------------------------------------------------------------//

    //The function space and some associate elements are then defined
    //Xh = space_type::New( mesh );

    element_struct_type U( M_Xh, "u" );
    element_struct_type V( M_Xh, "v" );


    //----------------------------------------------------------------------------------------------------------//

    //    F = vector_ptrtype( M_backend->newVector( Xh ) );
    //    D = sparse_matrix_ptrtype( M_backend->newMatrix( Xh, Xh ) );

    boost::timer time;
    time.restart();

    vector_ptrtype R( M_backend->newVector( M_Xh ) );
    sparse_matrix_ptrtype J( M_backend->newMatrix(M_Xh,M_Xh) );

    // set up the non linear solver
    M_backend->nlSolver()->residual = boost::bind( &StVenantKirchhoff<Dim, Order, Entity>::updateResidual,
                                                   boost::ref( *this ), _1, _2 );
    M_backend->nlSolver()->jacobian = boost::bind( &StVenantKirchhoff<Dim, Order, Entity>::updateJacobian,
                                                   boost::ref( *this ), _1, _2  );


    vector_ptrtype Uvec( M_backend->newVector( M_Xh ) );
    *Uvec = U;

    this->updateResidual( Uvec, R );
    this->updateJacobian( Uvec, J );

    M_backend->nlSolve( J, Uvec, R, 1e-10, 10 );
    U = *Uvec;

    cout << "[Solve] : " << time.elapsed() << "\n";

    this->exportResults( U, V );

    std::cout<< "\nNormL2 U :"<< U.l2Norm() << "\n";
  } // StVenantKirchhoff::run


} // Life


// instantiation
template class Life::StVenantKirchhoff<2,2,Life::Simplex>;
//template class Life::StVenantKirchhoff<3,2,Life::Simplex>;


